#include "../include/Distributions.hpp"
#include "../include/Density.hpp"
#include "../include/Kernel.hpp"
#include "../include/Particle.hpp"
#include "../include/Resample.hpp"
#include "../include/utils.hpp"
#include "../include/Models/Likelihood.hpp"
#include "../include/Misc/EP.hpp"
#include "../include/Misc/Laplace.hpp"
#include "../include/Misc/IWLS_Gelman.hpp"
#define NU 1 
extern "C"{

void IBIS(int *M, int *n, int *p, double *Y, double *pX,double *res, double *W,double *phi, double *acc)
{
	mat X=vect2mat<double>(pX,n[0],p[0]);
//	cout << X;
	//Density::Density_Probit *DP= new Density::Density_Probit(Y,X);
	mat m(p[0],1);
	mat foo(p[0],p[0]);
	m.zeros();
	foo.fill(0.1);
	mat s=diagmat(foo);
	
	mat bar(p[0],p[0]);
	bar.fill(2.5);
	mat S=diagmat(bar);

	//EP fEP;
	//fEP(X,Y,1/2.5,3);

//	Laplace fL;
//	fL(X,Y,m,S);	
	IWLS_G fIG;
	mat Sig;
	mat Mu;
	try{
		fIG(X,Y,m,S,NU);	
		Sig=fIG.Get_Sig();	
		Mu=fIG.Get_Mu();//mat Sig=S;
	}
	catch(std::exception& e)
	{
		cout << "Algorithm had numerical issue... Switching to simulating from prior";
		Sig=S;
		Mu=m;
	}
		//cout << fIG.Get_Mu();
	cout << Sig.diag();	
	//	cout << s;
	Distribution::Student *G;
	//Distribution::Student *G;
	//Distribution::Student *F=new Distribution::Student(p[0],NU,Mu,Sig);
	Distribution::Student *F;
	//Likelihoods::Probit3 *DP;
	Probit_st::Probit3 *DP;
	//Likelihoods::Probit2 *DP= new Likelihoods::Probit2(Y,X);
//	Kernel_MH<Density::IBIS,Distribution::Gaussian> *K= new Kernel_MH<Density::IBIS,Distribution::Gaussian>(DP,G);
	//Kernel_MH<Density::Tempering,Distribution::Gaussian> *K= new Kernel_MH<Density::Tempering,Distribution::Gaussian>(DP,G);
	Kernel_MH<Density::GeomBridge,Distribution::Student> *K;//	Kernel_MH<Density::GeomBridge,Distribution::Student> *K= new Kernel_MH<Density::GeomBridge,Distribution::Student>(DP,G);
	//Resample::SysResample *R=new Resample::SysResample;
	Resample::SysResample *R;
	//Particle<Kernel_MH<Density::IBIS,Distribution::Gaussian>,Resample::Multinomial,Density::IBIS> P(K,R,DP,M[0]);
	//Particle<Kernel_MH<Density::Tempering,Distribution::Gaussian>,Resample::SysResample,Density::Tempering> P(K,R,DP,M[0]);
//	Particle<Kernel_MH<Density::GeomBridge,Distribution::Student>,Resample::Multinomial,Density::GeomBridge> P(K,R,DP,M[0],F);
	Particle<Kernel_MH<Density::GeomBridge,Distribution::Student>,Resample::SysResample,Density::GeomBridge> *P;
//	cout << "hello";
	int numCPU = sysconf( _SC_NPROCESSORS_ONLN );
	numCPU=1;
	cout << "cpu  " << numCPU;
	mat A((int)(M[0]/numCPU),p[0]);
	mat w((int)(M[0]/numCPU),1);
	vector<mat> weights(numCPU,w);
	vector<mat> Pi(numCPU,A);
	vector<mat> Phi;
	vector<mat> Acc;
	mat Z(numCPU,1);	
	int i=0;
	//#pragma omp parallel for firstprivate(R,K,G,F,DP,P), shared(Pi,weights,Z, Phi,Acc)
		for(int i=0;i<numCPU;i++)
		{	
			//G=new Distribution::Gaussian(p[0],m,s);
			G=new Distribution::Student(p[0],NU,m,s);
			F=new Distribution::Student(p[0],NU,Mu,Sig);
			DP=new Probit_st::Probit3(Y,X,F,S,NU);
			//DP=new Likelihoods::Probit3(Y,X,F,S);
			K=new Kernel_MH<Density::GeomBridge,Distribution::Student>(DP,G);
			R=new Resample::SysResample;
			P=new Particle<Kernel_MH<Density::GeomBridge,Distribution::Student>,Resample::SysResample,Density::GeomBridge>(K,R,DP,(int)(M[0]/numCPU),F);
			P->Filter();
			Pi[i]=P->Get_theta();
			weights[i]=P->Get_W();
			Z(i,0)=P->Get_Z();	
			Phi.push_back(P->Get_Phiv());
			Acc.push_back(P->Get_Acc());
			delete F;
			//delete DP;
			delete G;
			delete R;
			delete P;
			delete K; 
		}
	//	cout << Phi[0];
	//	cout << Phi[1];
		mat Wi=vectortomat(weights,numCPU);
		mat theta=vectortomat(Pi,numCPU);
		
		mat Phim=vectortomat(Phi,numCPU);
		mat ACC=vectortomat(Acc,numCPU);

		//cout << Wi;	
//	cout << DP->Density(A.row(0).t(),0,n[0]-1);	
//	cout << A;
//	cout << "hello";
//cout << M[0] << p[0];	
	//double *w=P.Get_W();
//	*W=*w;
//
	cout << Z;
	double z=(as_scalar(mean(Z)));
	cout << "evidence: " << z;
	mat2vect<double>(theta,M[0],p[0],res);
	mat2vect<double>(Wi,M[0],1,W);
	mat2vect<double>(Phim,Phim.n_rows,1,phi);
	mat2vect<double>(ACC,ACC.n_rows,1,acc);

	/*	
	delete F;
//	delete DP;
	delete G;
	delete R;
	delete K;
	//delete P;	
*/
}

}
